###################################################################
## Extension Codes for Original Results
## Paper: Cities, Redistribution, and Authoritarian Regime Survival
###################################################################
rm(list=ls(all=TRUE))
library(readstata13)
library(stargazer)
library(MASS)
library(plm)
library(glmmML)
library(glmmADMB)
library(pglm)
library(survival)
library(ggplot2)
library(ggfortify)
library(sandwich)
library(lmtest)
library(mvtnorm)
library(foreach)
library(do)
library(maxLik)
library(ebal)
library(Matching)
library(optmatch)
setwd("/Users/1/Downloads/2018 Spring Harvard/Gov2001/Replication/2001_paper_docs")

# load the banks data
banks <- read.dta13("X-Ntl_Banks.dta")
banks$gdpgrow <- banks$gle_gdpgrowth/100
banks$ass <- banks$s17f1
banks$strike <- banks$s17f2
banks$gwar <- banks$s17f3
banks$govcri<- banks$s17f4
banks$purge <- banks$s17f5
banks$riot <- banks$s17f6
banks$rev <- banks$s17f7
banks$antidemon<-banks$s18f1
banks$surv.obj <- with(banks, Surv(time = `_t0`, time2 = `_t`, event = regimefail))
banks$external <- banks$strike + banks$gwar + banks$riot + banks$antidemon #+ banks$rev
banks$internal <- banks$ass + banks$govcri + banks$purge

banks$internal_0 <-banks$internal
banks$internal_0[banks$internal > 0 ] <- 1
banks$external_0 <-banks$external
banks$external_0[banks$external > 0 ] <- 1

# subsetting the datasets
military<- banks[which(banks$wr_mir ==1),]
singleparty <- banks[which(banks$wr_spr ==1),]
mornarchy <- banks[which(banks$wr_mor ==1),]
personal <- banks[which(banks$wr_pr ==1),]

######################################################################
######################################################################
######################################################################
######################## impacts of regime type on elite-level collective action
mod11 <- glm.nb(internal ~ (wr_mir + wr_mor + wr_spr)*(external)
                 + gdpgrow + gle_rgdp + urbcon +  r_atlas 
                 + conflict + leg + econcrisispwt
                 + factor(ccode) + factor(year),
                 data = banks)
summary(mod11)


mod12 <- glm(internal_0 ~ (wr_mir + wr_mor + wr_spr)*(external_0)
             + gdpgrow + gle_rgdp + urbcon +  r_atlas 
             + conflict + leg+ coldwar + econcrisispwt
             + factor(ccode) + factor(year),
              family=binomial(link='logit'), data=banks)
summary(mod12)

##################################### interpretations
####################### mod11 (nb)
set.seed(02138)
sim_combine <- rmvnorm(5000, mean = coef(summary(mod11))[,1],
                       sigma = vcovHC(mod11))
theta_sim <- rnorm(5000, mean = mod11$theta, sd = mod11$SE.theta)
sim_combine <- cbind(theta_sim, sim_combine)
dim(sim_combine)

X <- model.matrix(internal ~ (wr_mir + wr_mor + wr_spr)*(external)
                  + gdpgrow + gle_rgdp + urbcon +  r_atlas 
                  + conflict + leg+ coldwar + econcrisispwt
                  + factor(ccode) + factor(year), data = banks)

check <- as.vector(rownames(coef(summary(mod11))))
check
same <- colnames(X)[which(colnames(X) %in% check)]
X <- X[,same]
dim(X)
X <- na.omit(X)

median_cov <- rbind(apply(X, 2, median))
median_cov[1:10] # personalistic regime

result1 <- apply(sim_combine, 1, function(A){
  return(rnbinom(1, size = A[1],  mu = exp(median_cov%*%matrix(A[2:ncol(sim_combine)],ncol=1))))
})
mean(result1)
sd(result1)

median_cov2 <- median_cov 
median_cov2[2]<- 1 # military
result2 <- apply(sim_combine, 1, function(A){
  return(rnbinom(1, size = A[1],  mu = exp(median_cov2%*%matrix(A[2:ncol(sim_combine)],ncol=1))))
})
mean(result2)

median_cov3 <- median_cov
median_cov3[3]<- 1 #monar
result3 <- apply(sim_combine, 1, function(A){
  return(rnbinom(1, size = A[1],  mu = exp(median_cov3%*%matrix(A[2:ncol(sim_combine)],ncol=1))))
})
mean(result3)

median_cov4 <- median_cov
median_cov4[4]<- 1 # single-party
result4 <- apply(sim_combine, 1, function(A){
  return(rnbinom(1, size = A[1],  mu = exp(median_cov4%*%matrix(A[2:ncol(sim_combine)],ncol=1))))
})
mean(result4)

type <- rep(c("personalistic","military","monarchic","single party"), each = 5000)
table1 <- cbind.data.frame(c(result1, result2, result3, result4), type)
colnames(table1)[1]<- "result"
plot1 <- ggplot(table1,  aes(x = result, color = type)) + geom_density(adjust = 10, size = 0.8) + xlim(c(0,10))
plot1 + theme_bw()+labs(x = "frequency of elite-level collective action", y = "density")
mean(result2) - mean(result3)
mean(result2) - mean(result1)

## adding external collective action
df1 <- matrix(nrow = 11, ncol = length(median_cov))

for (i in 1:nrow(df1)){
  df1[i,] <- median_cov
}
colnames(df1) <- rownames(coef(summary(mod11)))
df1[,"external"] <- 0:10 #personalistic regime 
df2 <- df1 
df2[,"wr_spr"] <- rep(1, nrow(df2))
df2[,"wr_spr:external"] <- 0:10 #singe-party 
df <- rbind(df1, df2)

result <- c()
for (i in 1:nrow(df)){
  result1 <- apply(sim_combine, 1, function(A){
    return(rnbinom(1, size = A[1],  mu = exp(df[i,]%*%matrix(A[2:ncol(sim_combine)],ncol=1))))
  })
result <- rbind(result, result1)
}

diff <- apply(result,1,mean)
diff 
ex_value<- rep(0:10,2)
ex_value
type <- rep(c("personalistic","single-party"), each = 11)
type 
edge <- apply(result, 1, function(A){quantile(A, 1, probs = c(0.025, 0.975))})
edge <- t(edge)

table1 <- cbind.data.frame(diff,edge,ex_value, type)
colnames(table1)[2:3]<- c("lower","upper")
plot1 <- ggplot(table1) + geom_point(aes(x = ex_value,y = diff,color = type),pch = 18) + 
  geom_line(aes(x = ex_value,y = diff,color = type)) + 
  geom_ribbon(aes(x = ex_value, ymin=lower, ymax=upper, fill = type),alpha=0.3)
plot1 + theme_bw()+labs(x = "mass-level collective action", y = "elite-level collective action")

#### mod12
set.seed(02138)
sim_combine <- rmvnorm(5000, mean = coef(summary(mod12))[,1],
                       sigma = vcovHC(mod12))
dim(sim_combine)
X <- model.matrix(internal_0 ~ (wr_mir + wr_mor + wr_spr)*(external_0)
                  + gdpgrow + gle_rgdp + urbcon +  r_atlas 
                  + conflict + leg+ coldwar + econcrisispwt
                  + factor(ccode) + factor(year), data = banks)
check <- as.vector(rownames(coef(summary(mod12))))
same <- colnames(X)[which(colnames(X) %in% check)]
X <- X[,same]
dim(X)
X <- na.omit(X)

median_cov <- rbind(apply(X, 2, median))
median_cov[1:10] # personalistic regime

result1 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result1)
sd(result1)

median_cov2 <- median_cov 
median_cov2[2]<- 1 # military
result2 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov2%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result2)
sd(result2)

median_cov3 <- median_cov
median_cov3[3]<- 1 #monar
result3 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov3%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result3)
sd(result3)
mean(result3)

median_cov4 <- median_cov
median_cov4[4]<- 1 # single-party
result4 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov4%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result4)
sd(result4)
mean(result4)

prob1 <- apply(result1,2,mean)
prob2 <- apply(result2,2,mean)
prob3 <- apply(result3,2,mean)
prob4 <- apply(result4,2,mean)

type <- rep(c("personalistic","military","monarchic","single party"), each = 5000)
table2 <- cbind.data.frame(c(prob1, prob2, prob3, prob4), type)
colnames(table2)[1]<- "result"
plot2 <- ggplot(table2,  aes(x = result, color = type)) + geom_density(size = 0.8) + xlim(c(-0.01,1))
plot2 + theme_bw()+labs(x = "possibility of elite-level collective action", y = "density")
mean(result2) - mean(result3)
mean(result2) - mean(result1)

######################## impacts of regime type on mass-level collective action
#banks.new <- rbind.data.frame(military, singleparty, mornarchy, personal)
mod13 <- glm.nb(external ~  (wr_mir + wr_mor + wr_spr)*(internal)
                 + gdpgrow + gle_rgdp  + urbcon +  r_atlas
                 + conflict + leg+ coldwar + econcrisispwt
                 + factor(ccode) + factor(year),
                 data = banks)
summary(mod13)

mod14 <- glm(external_0 ~ (wr_mir + wr_mor + wr_spr)*(internal_0)
             + gdpgrow + gle_rgdp  + urbcon +  r_atlas
             + conflict + leg+ coldwar + econcrisispwt
             + factor(ccode) + factor(year),
             family=binomial(link='logit'), data=banks)
summary(mod14)




################ illustrations
set.seed(02138)
sim_combine <- rmvnorm(5000, mean = coef(summary(mod13))[,1],
                       sigma = vcovHC(mod11))
theta_sim <- rnorm(5000, mean = mod13$theta, sd = mod13$SE.theta)
sim_combine <- cbind(theta_sim, sim_combine)
dim(sim_combine)

X <- model.matrix(external ~ (wr_mir + wr_mor + wr_spr)*(internal)
                  + gdpgrow + gle_rgdp + urbcon +  r_atlas 
                  + conflict + leg+ coldwar + econcrisispwt
                  + factor(ccode) + factor(year), data = banks)

check <- as.vector(rownames(coef(summary(mod13))))
check
same <- colnames(X)[which(colnames(X) %in% check)]
X <- X[,same]
dim(X)
X <- na.omit(X)

median_cov <- rbind(apply(X, 2, median))
median_cov[1:10] # personalistic regime


## adding external collective action
df1 <- matrix(nrow = 11, ncol = length(median_cov))

for (i in 1:nrow(df1)){
  df1[i,] <- median_cov
}
colnames(df1) <- rownames(coef(summary(mod13)))
df1[,"internal"] <- 0:10
df2 <- df1 
df2[,"wr_mor"] <- rep(1, nrow(df2))
df2[,"wr_mor:internal"] <- 0:10 #mornar
df3 <- df1 
df3[,"wr_mir"] <- rep(1, nrow(df3))
df3[,"wr_mir:internal"] <- 0:10 #military 
df <- rbind(df1, df2, df3)
df4 <- df1 
df4[,"wr_spr"] <- rep(1, nrow(df4))
df4[,"wr_spr:internal"] <- 0:10 #military 
df <- rbind(df1, df2, df3, df4)

result <- c()
for (i in 1:nrow(df)){
  result1 <- apply(sim_combine, 1, function(A){
    return(rnbinom(1, size = A[1],  mu = exp(df[i,]%*%matrix(A[2:ncol(sim_combine)],ncol=1))))
  })
  result <- rbind(result, result1)
}
dim(result)

result.select <- c(result[6,],result[17,], result[28,], result[39,]) # when internal = 5
type <- rep(c("personalistic","monarchic","military","single party"), each = 5000)
table3 <- cbind.data.frame(result.select, type)
plot3 <- ggplot(table3,  aes(x = result.select, color = type)) + geom_density(adjust = 10, size = 0.8) + xlim(c(0,10))
plot3 + theme_bw()+labs(x = "frequency of mass-level collective action", y = "density")
mean(result[17,])-mean(result[28,])

################# glm
set.seed(02138)
sim_combine <- rmvnorm(5000, mean = coef(summary(mod14))[,1],
                       sigma = vcovHC(mod14))
dim(sim_combine)
X <- model.matrix(external_0 ~ (wr_mir + wr_mor + wr_spr)*(internal_0)
                  + gdpgrow + gle_rgdp + urbcon +  r_atlas 
                  + conflict + leg+ coldwar + econcrisispwt
                  + factor(ccode) + factor(year), data = banks)
check <- as.vector(rownames(coef(summary(mod14))))
X <- X[,check]
dim(X)
X <- na.omit(X)

median_cov <- rbind(apply(X, 2, median))
median_cov[1:10]# personalistic regime
colnames(median_cov)
median_cov[,"internal_0"] <- 1

result1 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result1)
sd(result1)

median_cov2 <- median_cov 
median_cov2[,"wr_mir"]<- 1 # military
median_cov2[,"wr_mir:internal_0"]<- 1 # military
result2 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov2%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result2)
sd(result2)

median_cov3 <- median_cov 
median_cov3[,"wr_mor"]<- 1 # military
median_cov3[,"wr_mor:internal_0"]<- 1 # military
result3 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov3%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result3)
sd(result3)
mean(result3)

median_cov3 <- median_cov 
median_cov3[,"wr_spr"]<- 1 # single-party
median_cov3[,"wr_spr:internal_0"]<- 1 # single-party
result4 <- apply(sim_combine, 1, function(A){
  odds = exp(median_cov4%*%matrix(A,ncol=1))      
  p    = odds / (1+odds)
  return(rbinom(n=4000, size=1, prob=p))
})
mean(result4)
sd(result4)
mean(result4)

prob1 <- apply(result1,2,mean)
prob2 <- apply(result2,2,mean)
prob3 <- apply(result3,2,mean)
prob4 <- apply(result4,2,mean)

type <- rep(c("personalistic","military","monarchic","single party"), each = 5000)
table4 <- cbind.data.frame(c(prob1, prob2, prob3, prob4), type)
colnames(table4)[1]<- "result"
plot4 <- ggplot(table4,  aes(x = result, color = type)) + geom_density(size = 0.8) + xlim(c(-0.01,1))
plot4 + theme_bw()+labs(x = "possibility of mass-level collective action", y = "")
mean(result2) - mean(result3)
mean(result3) - mean(result4)

stargazer(mod11,mod12, 
          title = "\\textbf{Regime Types Affect The Probability of Elite-Level Collective Actions}",
          keep = c("wr_mir",  "wr_mor", "wr_spr","external","external_0",
                   "gdpgrow", "gle_rgdp", "urbcon", "r_atlas", "conflict", 
                   "leg", "coldwar","econcrisispwt",
                   "wr_mir:external", "wr_mor:external", "wr_spr:external",
                   "wr_mir:external_0", "wr_mor:external_0", "wr_spr:external_0","Constant"),
          dep.var.caption = "",
          dep.var.labels = "",
          model.numbers = FALSE,
          no.space = TRUE,
          star.char = c("*", "**", "***", ""),
          column.labels = c("\\textbf{Negative Binomial}", "\\textbf{Logistic}"),
          covariate.labels = c("Military Regime", "Monarchic Regime", "Single-Party Regime",
                               "Mass-level Collective Action", "Mass-level Collective Action Dummy",
                               "GDP Growth (\\%)","Real GDP per Capita (logged)",
                               "Urban Concentration","Ethnolinguistic Fractionalization",
                               "Civil/Intl. Conflict","Legislature", "Cold War","Economic Crisis",
                               "Military Regime*Mass-level Collective Action",
                               "Monarchic Regime*Mass-level Collective Action", 
                               "Single-Party Regime*Mass-level Collective Action",
                               "Military Regime*Mass-level Collective Action Dummy",
                               "Monarchic Regime*Mass-level Collective Action Dummy", 
                               "Single-Party Regime*Mass-level Collective Action Dummy","Constant"),
          omit.stat = c("rsq", "aic", "ll", "theta"),
          notes = "\\parbox[t]{\\textwidth}{\\textit{Note:} 
          *** p $<$ 0.01, ** p $<$ 0.05, * p $<$ 0.1.}",
          notes.align = "l",
          notes.append = FALSE,
          notes.label = "")


stargazer(mod13,mod14, 
          title = "\\textbf{Regime Types Affect The Probability of Mass-Level Collective Actions}",
          keep = c("wr_mir",  "wr_mor", "wr_spr","internal","internal_0",
                   "gdpgrow", "gle_rgdp", "urbcon", "r_atlas", "conflict", 
                   "leg", "coldwar","econcrisispwt",
                   "wr_mir:internal", "wr_mor:internal", "wr_spr:internal",
                   "wr_mir:internal_0", "wr_mor:internal_0", "wr_spr:internal_0","Constant"),
          dep.var.caption = "",
          dep.var.labels = "",
          model.numbers = FALSE,
          no.space = TRUE,
          star.char = c("*", "**", "***", ""),
          column.labels = c("\\textbf{Negative Binomial}", "\\textbf{Logistic}"),
          covariate.labels = c("Military Regime", "Monarchic Regime", "Single-Party Regime",
                               "Elite-level Collective Action", "Elite-level Collective Action Dummy",
                               "GDP Growth (\\%)","Real GDP per Capita (logged)",
                               "Urban Concentration","Ethnolinguistic Fractionalization",
                               "Civil/Intl. Conflict","Legislature", "Cold War","Economic Crisis",
                               "Military Regime*Elite-level Collective Action",
                               "Monarchic Regime*Elite-level Collective Action", 
                               "Single-Party Regime*Elite-level Collective Action",
                               "Military Regime*Elite-level Collective Action Dummy",
                               "Monarchic Regime*Elite-level Collective Action Dummy", 
                               "Single-Party Regime*Elite-level Collective Action Dummy","Constant"),
          omit.stat = c("rsq", "aic", "ll", "theta"),
          notes = "\\parbox[t]{\\textwidth}{\\textit{Note:} 
          *** p $<$ 0.01, ** p $<$ 0.05, * p $<$ 0.1.}",
          notes.align = "l",
          notes.append = FALSE,
          notes.label = "")

######################## 
######################## impacts of collective action on regime survival
######## A general look 
mod <- coxph(surv.obj ~  ass + strike + gwar + govcri + purge + riot + rev + antidemon + r_atlas + econcrisispwt + coldwar  + urbcon + leg + conflict + gdpgrow + gle_rgdp,
             data = banks, ties = "efron",
             robust = TRUE)
summary(mod)
est <- mod$coefficients[1:8]
exp(est)
se <- coef(summary(mod))[1:8,4]

y.name <- c("assassinations","general strikes","guerilla warfare","government crises",
            "purges","riots","revolutions","anti-gov demonstrations")

table1 <- cbind.data.frame(y.name,est,se)
plot1 <- ggplot(table1,aes(y = y.name, x = est)) + geom_point(pch=18) + 
  geom_errorbarh(aes(xmin=est-1.96*se, xmax=est+1.96*se),height = 0.4, color = "grey40") +
  geom_vline(xintercept =0,col = "indianred",lty = 2)
plot1+ theme_bw()+labs(x = "estimates", y = "collective action")

modmod <- coxph(surv.obj ~  ass + strike + gwar + govcri + purge + riot + rev + antidemon + gdpgrow + gle_rgdp,
                data = banks, ties = "efron",
                robust = TRUE) ## trying different covariates
summary(modmod)

######################## different regimes 
# collective action to durability
mod21 <- coxph(surv.obj ~  internal + external + r_atlas + econcrisispwt + coldwar  + urbcon + leg + conflict + gdpgrow + gle_rgdp,
               data = military, ties = "efron",
               robust = TRUE)
summary(mod21)

mod22 <- coxph(surv.obj ~ internal + external  + r_atlas + econcrisispwt + coldwar + urbcon + leg + conflict + gdpgrow + gle_rgdp,
               data = singleparty, ties = "efron",
               robust = TRUE)
summary(mod22)

mod23 <- coxph(surv.obj ~ internal + external + r_atlas + econcrisispwt + coldwar + urbcon + leg  + conflict + gdpgrow + gle_rgdp,
               data = mornarchy, ties = "efron",
               robust = TRUE)
summary(mod23)

mod24 <- coxph(surv.obj ~ internal + external  + r_atlas + econcrisispwt + coldwar + urbcon + leg  + conflict + gdpgrow + gle_rgdp,
               data = personal[personal$gle_pop > 1000, ], ties = "efron",
               robust = TRUE)
summary(mod24)

robust.se <- list(sqrt(diag(mod21$var)),
                  sqrt(diag(mod22$var)),
                  sqrt(diag(mod23$var)),
                  sqrt(diag(mod24$var)))

stargazer(mod21, mod22, mod23, mod24, se = robust.se,
          title = "\\textbf{Mass-/Elite- Level Collective Action Harms Regime Survival}",
          dep.var.caption = "",
          dep.var.labels = "",
          model.numbers = FALSE,
          no.space = TRUE,
          star.char = c("*", "**", "***", ""),
          column.labels = c("\\textbf{Military}", "\\textbf{Single-Party}",
                            "\\textbf{Monarchic}","\\textbf{Personalist}"),
          covariate.labels = c("Elite-level Collective Action", "Mass-level Collective Action",
                               "Ethnolinguistic Fractionalization","Economic Crisis","Cold War",
                               "Urban Concentration", "Legislature",
                               "Civil/Intl. Conflict","GDP Growth (\\%)",
                               "Real GDP per Capita (logged)"),
          omit.stat = c("rsq", "max.rsq", "wald", "lr", "logrank"),
          #add.lines = list(c("Failures", mod21$nevent, mod22$nevent,mod23$nevent)),
          notes = "\\parbox[t]{\\textwidth}{\\textit{Note:} 
          All models in this table are Cox Models. Positive
          coefficients reflect estimated effects that increase
          hazard rates of collapse. Robust standard 
          errors in parentheses. *** p $<$ 0.01, ** p $<$ 0.05, * p $<$ 0.1.}",
          notes.align = "l",
          notes.append = FALSE,
          notes.label = "")

########################################
########################################
########################################
########################################

